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Abstract 

The steady-state and nonequilibrium properties of the model of environmental- 
economic interactions are studied. The interacting heterogeneous agents are simu- 
lated on the platform of the emission dynamics of cellular automaton. The model 
possess the discontinuous transition between the safe and catastrophic ecology. 
Right at the critical line, the broad-scale power-law distributions of emission rates 
have been identified. Their relationship to Zipf's law and models of self-organized 
criticality is discussed. 
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1 Introduction 



On the onset of global ecological crisis questions related to regulation technolo- 
gies and safe handling of environmental resources become urgent. The equilib- 
rium between the economic growth and environmental variables is an actual 
problem of every government. The environmental science must deal with the 
task how to design robust sensorial nets and control systems to manage the 
balance between the immediate economic needs and the long-time perspectives 
of mankind. There are several instruments for control and decrease the local 
and global pollution rate: energy taxation, pollution rights and introduction 
of cleaner technologies. Unfortunately, direct and naive application of these 
approaches does not result relevant emission reductions. Enforcing pollution 
rights is difficult and has high administrative costs [1,2]. As an example can 
serve the multi-country control model of international coordination of emis- 
sion charges and pollution trading [3]. As shown in Ref. [3] pollution trading 
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may result a non-cooperative equilibrium with high pollution rates. It is often 
argued that if the only tool of policy is the energy taxation it may cause the 
loss of the competitiveness between firms. Moreover, the tax rate which is be 
able to stabilize the emissions would be too high to be economically and po- 
litically acceptable. Because of energy taxation and network effects in society, 
the firms tend to delay adoption to cleaner technologies [4,5]. All these unex- 
pected emergent regulatory effects stimulated our interest about the problem 
of dynamics of emissions that can be seen as a consequence of the competitive 
interaction of the environmental and economical (acquisition) agents. 

A key principle we applied to define a robust environmental policy is the 
feedback mechanism. In general, the nonlinear feedback is one of the key up- 
bottom mechanism that enables to keep control over the dynamics of complex 
systems. The mechanism has been exploited in many different systems. The 
problem discussed in Ref . [6] is focused to the feedback localization of the para- 
metric boundary between periodicity and chaos generated by a logistic map. 
Stabilization of dynamical regimes via feedback appears in the simulations of 
investments Ref. [7]. In Ref. [8,9] model has been suggested, where the specula- 
tive market expectations are linked via the positive feedback to the prevailing 
price trend. As stated in Ref. [10], the steady self- organized critical (SOC) 
regime Ref. [11] implies the operation of inherent feedback mechanism that 
ensures a steady state marginally stable against disturbances. It is important 
here to note that criticality is accompanied by the absence of characteristic 
scale of the power-law distributions of avalanches. According Ref. [12] the non- 
linear feedback allows to transform the models of the "unstable" phase transi- 
tions (i.e. classic models) to SOC models. In the works devoted to self-adjusted 
Monte Carlo methods Ref. [13,14], the feedback which drives the physical sys- 
tem to a critical regime is expressed in terms of estimated statistical averages. 
Additionally, catastrophic environmental collapses in consumption/pollution 
feedback models were described in [15,16]. In this paper we present feedback 
control distributed over the population of autonomous sensing mobile agents 
that control the behavior of factory agents. 

As usual, the thinking about distributed control is motivated by the expected 
robustness against accidental failures. Nevertheless, most of agent systems 
face to the general problem of agency that is the imperfect incentives [17] to 
satisfy the global goal function of system. Our motivation is to show, that this 
property of distributed agent system may drive the environment to a critical 
state with extremal emission ratio. The idea of walking sensorial agents in our 
model is very similar to the idea of distributed pollution sensing by pigeon 
bloggers proposed by de Costa Ref. [18], where pigeon agents are equipped by 
communication chips and sensors for carbon monoxide. Similar idea had been 
reported as a concept of randomly walking sensors for detection of dynamically 
changing gradient sources [19]. 
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The plan of the paper is as it follows. In the next section the multi-agent 
ecological model is introduced. In section sec. 3 the results of numerical study 
are discussed. Finally the conclusions are presented. 



2 The model 

The first point we consider is the dynamics of diffusive emissions of factory 
agents. The base substance of model is a square lattice with dimension L x L 
which represents the physical space for emission spreading. The factories are 
located on random lattice sites and the spatial distribution of their emissions is 
directed by discrete rule of diffusion. For this purpose, we implement a cellular 
automaton (CA) based on integer rules. Our suggestion for computation of 
diffusion is partially inspired by the lattice gas models of pollutants Ref.[20]. 
However, we preferred formulation where the concentration field on lattice 
sites varies via the recursive rule 

VRG{LxL}, m(f + l,R) =m(t,R) + / nn (t,R) + / F (t,R) (1) 

that replaces in one time step t the m(t, R) "integer pieces" of emission at the 
position R = [R x , R y \. The emission inflow caused by the factory agents I F is 
discussed in below Eq.(6) The emission current between the nearest neighbor 
cells 

I nn (t, R) = -A[m(t, R)/5] int + J2 i m (^ r )/ 5 W , (2) 

renn(R) 

where the brackets [. . .] mt denote the integer part of argument and Z^nn(R) 
refers to the four nearest neighboring sites of R. Eq.(2) preserves the total 
emission for transitions between the lattice sites . Other words, the natural 
purification of emission is absent, emissions are exclusively depleted by the 
outflow through opened lattice bounds. For m(t, R) ^> 1 the mass equidis- 
tribution of emissions becomes satisfied with \j^pm order. Under these con- 
ditions, I nn (t, R) acts as 2d Laplace of m, and thus CA rules describe the 
stylized diffusive transport. 

The local concentration of emissions is permanently monitored by N§ number 
of sensing agents. The sensing agents are random walkers with the lattice 
positions S(t,j) = [S x (t,j), S y (t,j)], j = 1,2, . . ., N$. The lattice position of 
the jth walker is given by 

S x (t + l,j) = [S x (t,j) + s x (t,j)] modL , (3) 
S y(t + 1 J) = l S y(tJ) + s y(tJ)} mod L j 
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where the shift s(t,j) = [s x (t,j),s y (t,j)} is chosen randomly from the set of 
four unit vectors {(0, 1), (0, —1), (1, 0), (—1, 0)}. Eq.(3) preserves the periodic 
boundary conditions for S(t,j). Moreover, we checked that the statistical de- 
pendencies studied in this paper do not change remarkably for closed boundary 
conditions. 

The source of emissions are the factory agents localized at random positions 
F(t, k) e {L x L}, k = 1,2, ... , Np . They are characterized by an integer 
value of emission production n(t, k). The condition of sensorial feedback sys- 
tem is that over a predetermined local pollution threshold m c every sensorial 
agent with m(t, S(t,j)) > m c sends a feedback request to the nearest (in the 
Euclidean sense) factory F(t, /^iD 

min \\S(t,j)-F(t,k)\\ = \\S(t,j)-F(t,k£l)\\. (4) 

k=l,2,...,N F 



The concentration threshold m c is analogous to that in Ref.[21] where the 
feedback operation is mediated by the chemical sensors with specific thresholds 
for different pollutants. Any request of sensorial agent instantly decreases the 
local emission production n(t, k) of the A;th factory by a constant additive 
factor A > 0. In the absence of feedback signal, to satisfy economic needs, the 
factory increases its production automatically. Both alternatives are described 
by the equation 



Vfc, n(t + l,k) = < 



n(t,k) + A 



max{n(t, k) — A; 0} 



when no signal 

from any S(t, j) agent 

is received 

feedback signal 
received 



(5) 



The eventuality n(t, k) — A < (zero emission means no production) yields the 
bankruptcy of F(t, k) and the appearance of a new factory at random position 
with initial production n(t + 1, k) — 0. The zero initial emission assure that 
new factory may start its production only when it survives at least one time 
period without feedback request. The total number of factories is conserved 
by the death-birth rules. Since the production strategies of new generation 
of factories does not change, the implicit assumption here is that system is 
locked to the inferior technology [22]. 
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What remains unspecified in Eq.(l) is the emission inflow caused by the F 
type of agents 

I F (t, R) = n(t, k) 5 R ,F(t,k) ■ (6) 

It is worth nothing that all the agent systems in common with CA subsystem 
are updated synchronously. The asynchronous update rules can be also of 
interest, but we dismiss them from our considerations. 

After the completion of the agent-based model, its relationship to the stan- 
dard sandpile model of SOC [11] should be noted. The level of emissions is 
confronted with a threshold which is similar to the height (slope) threshold 
in SOC models. In sandpile models if the height of a pile exceeds a given 
threshold, a "signal" is transmitted to the central site and to its neighboring 
sites. The abundant grains are then toppled from the center. Spatial expan- 
sion of a sandpile avalanche corresponds to the emission spreading around a 
factory agent. The power-law distribution of avalanches in sandpile model is 
a consequence of the local conservation of sand grains. The property coincide 
with the conservation of emissions, but contradicts to the non-local conserva- 
tion of factory agents. The opened boundary conditions for emissions are also 
borrowed from the sandpile model. 



3 Simulation results 

The simulations have been performed for lattices L = 20, 30, 40, 50. To avoid 
the influence of transients, the initial 10 5 updates per site were discarded; 
about r = 5 x 10 7 updates have been used to collect information about the 
steady-state for the fixed parameters m c — 5 >> A = 1. In further we limit 
ourselves to the systems with p-p = Np/L 2 = 1/100. In the world with large 
number of sensors iVg the emissions fluctuate around m c . This oversaturated 
population of the sensors prevents from unbearable contamination but yields 
a non-efficient economy. Much different is the situation for N s small, where 
the steady-state regime become unstable since the rate of emissions diverges 
in time. The dynamics corresponds to a catastrophic regime. The preliminary 
numerical work on steady-state statistics yields the conclusion that m(t, R) is 
remarkably lowered at the corners of L x L lattice world with opened bounds. 
Such inhomogeneity induces secondary effect of the accumulation of factories 
at boundary regions. As L increases, the overall contribution of corners decays 
faster than the contribution of edges. In this context, the important is the 
aspect of finite-size scaling carried out in the next section. 

The outcome of our simulation experience is the decision to accumulate the 
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statistics of those factories that have received even one feedback request during 
the step t. Thus, we define the following order parameter 

i t-1 i N F 

(0) = -£0(i), <Kt) = — Y,°(t,k), (7) 

T t=0 7V F fc=l 

where (0) is the temporal average calculated over r epochs. If source F(t, k) re- 
ceived a feedback signal a(t, k) — 1 and cr(£, fc) = otherwise. Fig.l shows the 
N$ dependencies of (0) that attain abrupt jump interpreted as a discontinuous 
transition form the safe (nondestructive) to the destructive economy. More- 
over, we found that N s dependencies of (0) collapse onto the single finite-size 
scaling form 

(0) = £-/(£), e = N s /L 2+b (8) 



with exponents a = 0.38, b = 1.39. The algebraic combination crit = N s /L 2+b , 
where the scaling function /(...) attains the jump /(0 CT it) an d (0) tends to 
the L~ a asymptotics, is interpreted as a critical point. For given pp we have 
estimated 6 cv \ t — 2.8 x 1CT 4 . For sufficiently small lattices (L < 20) the scaling 
violation is indicated. 

The problem with sensorial system manifests itself when 9 < 9 cr \ t . In this 
regime a sensor multiply affects the same nearest factory and lowers its emis- 
sion production, while other hidden sources, maybe more influential and dan- 
gerous polluters, could left without feedback requests. Hereafter, this process is 
referred as spatial screening. The screening emerges as a result of co-evolution 
and self- organization in multi-agent system. The statistical significance of the 
screening is clarified in Fig. 2 where the auxiliary order parameter 




is plotted against 9. The cr'(t,j) = 1 when the request comes from the jth 
sensorial agent, a'(t,j) = otherwise. At the critical point 9 = 9 ciit both 
the order parameters attain discontinuous jump. The ratio N^/Np reflects the 
screening efficiency which at the critical point reads N^/Np = 9 cr - lt L l '/ 1 p-p = 
9 C ritPp <yh+2 ^ 2 Np 2 . The proportionality N$/Np ~ Np 2 reveals the interpreta- 
tion of the index b. The difference (0') — (0) reflects the impact of the screening 
as well as the most pronounced situation at 9 < 9 CT - lt . 

The vicinity of criticality is corroborated by the broad-scale statistics of pro- 
duction rates (see Fig. 3). We see that the most pronounced power- law depen- 
dence belongs to 9 ~ 9 crit . In the economic models Ref.[3], the emission rates 



6 



0.22 
0.2 
0.18 
0.16 
0.14 
0.12 
0.1 
0.08 
0.06 
0.04 
0.02 



f- 




a = 0.38 



6= 1.39 



L 
L 
L 
L 



20 
30 
40 
50 







0.0001 0.0002 0.0003 0.0004 0.0005 



L 2+b 



Fig. 1. Finite size scaling of the order parameter (0) as a function of sensorial agent 
density Ns/L 2 scaled by L~ b factor. The simulation was performed for different L 
with fixed density of factories pp = 0.01. 

are usually related to the overall production of goods with a proportionality 
factor called emission-output ratio. Standard and widely accepted assumption 
is that this ratio is fixed and irrespective to the investments in clean technol- 
ogy. With such idealization, the critical fluctuations of n(t, k) can be assigned 
to the empirical Zipf 's law. As we see from Fig. 4, the cumulative probability 
distribution function pdf > (n) corresponding to (iV s ,L)- parametric domain 
have been characterized by the effective exponents within the broad range 
P ~ 0.5 — 6.8. It means that feedback mechanism we proposed is not sufficient 
to fix the realistic exponents of the Zipf's law (according Ref. [23] the empirical 
exponents are (5 = 0.84; 0.995). 

Since the place for new factory is chosen randomly, the snapshots showing the 
spatial distributions of F agents evoke the images of the short-range ordered 
molecular patterns. From this perspective the system of factories have fea- 
tures of 2d molecular liquid or "vapor". Within such viewpoint the "vapor" 
of F(t, k) agents prepared for 6 > 9 crit belongs to the safe economy, whereas 
stronger correlated "liquid condensate" (9 < 6 ciit ) corresponds to the catas- 
trophic regime. In further, these arguments are applied in favor of scenario 
of discontinuous transition. The transition concept is examined in further 
(in a manner that is independent of (0)) via damage spreading (DS) method 
Ref. [24] . The DS is based on the comparative analysis of simultaneously evolv- 
ing replicas. Let us consider the pair of the systems I, II with equal number of 
agents. At certain time t the copy II is replicated from I and, consequently, 
single grain of the emission (damage) is added to its randomly chosen cell R': 



VR, t' = 0, 



m 



( n )(0,R) 



= m 



(t, R) + <5 r ,r, . 



(10) 
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Fig. 2. The order parameters as function of sensorial agent density N$/L 2 scaled 
by L~ b factor for L = 40. The other parameters as in Fig.l. 
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Fig. 3. The log-log plot of the cumulative distribution of emissions n(t, k) for L = 50. 
The boundary power-law regime (iVs = 160) is the most pronounced. Thereupon 
for iVs < 160 the dynamics has no steady state and steady-state distributions do 
not exist. 



The subsequent t' time steps are executed with the identical sequences of 
random numbers for updating of the both of replicas I and II. The whole 
process starting from the configuration given by Eq.(10) is repeated many 
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Fig. 4. The cumulative pdf 's of the production rates n(t, k) obtained for different 
L and Ng(L) = 9 CT i t L b+2 The inset table lists the effective exponents (3. 
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Fig. 5. The application of DS technique for L = 40, t' = 250. The averaged 
inter-replica distance D(t') [see Eq.(ll)] is plotted versus 6 = N s /L 2+b . The position 
of the critical point coincides with Fig.l within the statistical error. 

times. The idea of DS is to look for averages (D(t')) of the Hamming distance 

(11) 



■y LxL 

D (t') = -j~2 (1 ~ ^m( I )(t+t',R),m( II )(t',R)) > 



R 



where D(0) = 1/L 2 . In Fig. 5 we show (D(t')) calculated as an average taken 
over the 1000 steps of the steady-state. It is seen that DS technique confirms 
^-driven transition 6 = 6 cr \t maintained previously by the direct inspection of 
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Fig. 6. The differences between nonequilibrium relaxation below and above the crit- 
ical number (i.e. minimum number) of sensors sufficient to prevent from dangerous 
pollution is iV s ^ 80 (for L = 40). The interval (2000 < t < 4000) approximately 
delimits the fast CA purification through the opened bounds. The steady-state be- 
comes formed above t ~ 7000 (above the crossing point). 

(</>)(#) dependence. 

Next, the system is studied via nonequilibrium relaxation method [25] that 
is an efficient tool of analysis of equilibrium phase transitions. The numerical 
results are shown in Fig. 6. The initial configuration has been prepared in a such 
way that at the beginning, the factories are generated at random positions. 
Subsequently, the emissions have been left to increase freely without control. 
After the level J2r xL r ) = 5000L 2 is reached, the sensors become active. 
The resulting relaxation dependence is then constructed by averaging over 
1000 repetitions of the whole relaxation process. We found that relaxation 
strongly depends on how far 6 is from 6* cr i t . Three distinct dynamical regimes 
have been identified when inspecting averaged relaxation process. In early 
stage of development practically all sensors are active since ((f)) ~ 1. Later, 
when content of the most of cells are still above m c , relatively fast exponential 
relaxation dominates. It corresponds to the purification through bounds. As 
it can be clearly seen, the processes at large times are well separated from 
the initial relaxation. The long-time process is associated with the positional 
adjustment (or self-organization) of the factories. 



4 Conclusions 

The statistical properties of the multi-agent system including the environmen- 
tal-economic interactions have been studied. It should be noted that at early 
stages of the work we were motivated by the particular hypothesis about im- 
plementation of feedback and direct generation of SOC regime. However, even 



10 



preliminary runs have uncovered that not all the regimes conducted by feed- 
back acquire SOC signatures as are the power-law distributions and power-law 
autocorrelation functions. With these facts in mind, we shift our focus to give 
an interpretation in terms of phase transitions ("unstable" in sense of Ref. [12]). 
Our simulation confirmed that only optimized feedback induces attractiveness 
of the critical state. In other words, the distributed feedback (as many other 
approaches) does not automatically guarantees the attractiveness of critical 
point. The agents proposed here hold only rudimentary thinking and thus more 
intelligent entities should be considered in future work. A variety of problems 
can be formulated where the establishment of new factory is subjected to a 
place-search optimization strategy. With such modification the formation of 
highly correlated factory - sensorial agent structures are expected. One can 
extend the idea to models with crystal-like ordering and formation of the 
compact aggregates of factories - industrial complexes enhancing the effect of 
screening. An alternative representation of distributed feedback mechanism 
may be the notion of individual energy taxation of factories which is pro- 
portional to the local contamination level. Contrary to plethora of possible 
additional realistic entries omitted, the hypothesis posed by simulation is that 
the phase transition appears to be a generic statistical aspect of the groups of 
autonomous agents with conflicting interests. The main implication from the 
present model is the identification of power-law (nearly catastrophic) pdf's of 
fluctuations of emission rates. We believe that focus to power-law distributions 
will open new perspectives of the environmental diagnostics. 
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